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Abstract 

A physical-mathematical approach to anomalous diffusion may be based 
on generalized diffusion equations (containing derivatives of fractional order 
in space or/and time) and related random walk models. By space-time 
fractional diffusion equation we mean an evolution equation obtained from 
the standard linear diffusion equation by replacing the second-order space 
derivative with a Riesz-Feller derivative of order a G (0,2] and skewness 
9 {\9\ < min{a,2 — a}), and the first-order time derivative with a 
Caputo derivative of order /? G (0, 1] . Such evolution equation implies for 
the flux a fractional Pick's law which accounts for spatial and temporal 
non-locality. The fundamental solution (for the Cauchy problem) of the 
fractional diffusion equation can be interpreted as a probability density 
evolving in time of a peculiar self-similar stochastic process that we view 
as a generalized diffusion process. By adopting appropriate finite-difference 
schemes of solution, we generate models of random walk discrete in space 
and time suitable for simulating random variables whose spatial probability 
density evolves in time according to this fractional diffusion equation. 
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1 Introduction 

It is well known that the fundamental solution (or Green function) for the 
Cauchy problem of the linear diffusion equation can be interpreted as a 
Gaussian (normal) probability density function (pdf) in space, evolving in 
time. All the moments of this pdf are finite; in particular, its variance 
is proportional to the first power of time, a noteworthy property of the 
standard diffusion that can be understood by means of an unbiased random 
walk model for the Brownian motion. 

In recent years a number of master equations have been proposed for 
random walk models that turn out to be beyond the classical Brownian 
motion, see e.g. Klafter et al. [32j. In particular, evolution equations 
containing fractional derivatives have gained revived interest in that they are 
expected to provide suitable mathematical models for describing phenomena 
of anomalous diffusion and transport dynamics in complex systems, see e.g. 
[U, [To], [E], [28] [Ml [3D], [HllMllaS], [40], [481149], [52], [57], [SH [62], 
[631 I64j . For a recent review we refer to Metzler and Klafter [H] where 
other references are found. 

Here we intend to present our original approach to the topic that, being not 
considered in [41], could offer some novel and inspiring inspections to the 
phenomenon of anomalous diffusion which is of great interest in chemical 
physics. In this paper we complement and revisit some of our previous 
results found e.g. in [19], [23], [26]. We first show that our proposed 
fractional diffusion equations can be derived from generalized Pick's laws 
which account for spatial and/or temporal non- locality. Then we pay 
attention to the fact that the fundamental solutions (or Green functions) of 
our diffusion equations provide spatial probability densities evolving in time, 
related to self-similar stochastic processes, that we view as generalized (or 
fractional) diffusion processes to be properly understood through suitable 
random walk models. More precisely, we replace the second-order space 
derivative or /and the first-order time derivative by a suitable integro- 
differential operator, which can be interpreted as a space or time derivative 
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of fractional order a G (0, 2] or /5 G (0, 1] , respectivelj0. The space 
fractional derivative is required to depend also on a real parameter d (the 
skewness) subjected to the restriction |^| < min{a, 2 — a} . Correspondingly, 
the generalized equation will be referred to as the strictly space fractional 
diffusion equation of order a and skewness ^ if a G (0, 2) and /? = 1 , or the 
strictly time fractional diffusion equation of order (3 if a = 2 and (3 G (0, 1). 
In general, allowing a G (0, 2) and [3 G (0, 1) , we have the strictly space-time 
fractional diffusion equation of order a, /3 and skewness 6 . Of course, in the 
case {a = 2, /5 = 1} we recover the standard diffusion which leads to the 
Gaussian probability density and to the classical Brownian motion. 

For the strictly space fractional diffusion of order a ({0 < a < 2, /3 = 1}) 
we generate the class of (non-Gaussian) Levy stable densities of index a and 
skewness 9 {\6\ < min{a, 2 — a}), according to the Feller parameterization. 
As known, these densities exhibit fat tails with an algebraic decay oc 
1^1 ^ thus obtain a special class of Markovian stochastic processes, 
called stable Levy motions, which exhibit infinite variance associated to the 
possibility of arbitrarily large jumps {Levy flights). 

For the strictly time fractional diffusion of order f3 {{a = 2, < /3 < 1}) 
we generate a class of symmetric densities whose moments of order 2n are 
proportional to the n/3 power of time. We thus obtain a class of non- 
Markovian stochastic processes (they possesses a memory!) which exhibit a 
variance consistent with slow anomalous diffusion. 

For the strictly space-time fractional diffusion of (composite) order a G 
(0, 2), /? G (0, 1) we generate a class of densities (symmetric or not symmetric 
according to 9 = or 6 0) which exhibit fat tails with an algebraic decay 
cx . Thus they belong to the domain of attraction of the Levy 

stable densities of index a and can be called fractional stable densities. The 
related stochastic processes, by possessing the characteristics of the previous 
two classes, are non-Markovian and exhibit infinite variance; however, 
the possibility of arbitrary large jumps is contrasted by memory effects. 
Furthermore we mention the cases a = f3 for which it is possible to derive 
the Green function in closed analytical form: we refer to these cases as to 
neutral diffusion. 



^We remind that the term "fractional" is a misnomer since the order can be a real 
number and thus is not restricted to be rational. The term is kept only for historical 
reasons, see e.g. ^2]. Our fractional derivatives are required to coincide with the standard 
derivatives of integer order as soon as a = 2 (not as a = 1!) and (3 — 1 . 
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We shall prove that in any case the corresponding Green function can be 
interpreted as a spatial probability density evolving in time with a self- 
similarity property having scaling exponent v = fija . This allows us to limit 
ourselves to consider the expression of the Green function at a fixed time, say 
t = 1 , namely to the so-called reduced Green function. To approximate the 
time evolution of all the above densities we propose finite difference schemes, 
discrete in space and time, for the fractional derivatives. By taking care in 
constructing these schemes, namely by requiring them to be conservative 
and non-negativity preserving, they can be interpreted as discrete random 
walk models for simulating particle paths by the Monte Carlo technique. By 
properly scaled transition to vanishing space and time steps, these models 
can be shown to converge to the corresponding continuous processed. 

The paper is divided as follows. In Section 2 we first present our space 
and time fractional diffusion equations providing the definitions of the 
space and time fractional derivatives based on their Fourier and Laplace 
representations, respectively. Then we show how to derive them from 
generalized Pick's laws. Section 3 is devoted to the Green functions, pointing 
out their similarity properties. We provide the representations of the 
corresponding reduced Green functions in terms of Mellin-Barnes integrals 
which allow us to obtain their computational expressions. In Section 4, we 
first discuss the discrete random walk approach to the Brownian motion, 
which is based on the well-known discretization of the second order space 
derivative and the first-order time derivative entering the standard diffusion 
equation. Then, by properly discretizing the space-fractional derivative we 
generalize the above approach to the more general Markovian case of strictly 
space fractional diffusion. Section 5 is devoted to the extension of the above 
approach to the non-Markovian cases of strictly time fractional diffusion 
and strictly space-time fractional diffusion equations. Section 6 is devoted 
to the numerical results of our random walks produced in some case-studies 
and to the concluding discussions. For possible convenience of the reader 
we have reserved Appendix A and Appendix B for treating with some detail 
the Riesz-Feller and Caputo fractional derivatives, respectively. 

^ This was shown by Gorenflo and Mainardi for the space fractional diffusion in [231 
1241 125] . For the general case it will be shown in a next paper. 
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2 The space-time fractional diffusion equation 



By replacing in the standard diffusion equation 
d 52 

—u{x,t) = —-^u{x,t), —oo<x<+(X), t>0, (2.1) 
ot ax^ 

where u = u{x, t) is the (real) field variable, the second-order space 
derivative and the first-order time derivative by suitable integro- differential 
operators, which can be interpreted as a space and time derivative of 
fractional order we obtain a generalized diffusion equation which may be 
referred to as the space-time fractional diffusion equation. We write this 
equation as 

tD^u{x,t) = a:DgU{x,t), -OO < X < +00 , t>0, (2.2) 

where a , 6 , fi are real parameters restricted to 

0<a<2, \e\ <mm{a,2 - a} , < /? < 1 . (2.3) 

In (2.2) is the Riesz-Feller fractional derivative (in space) of order a 

and skewness 6 , and is the Caputo fractional derivative (in time) of 
order /3 . 

The definitions of these fractional derivatives are more easily understood if 
given in terms of Fourier and Laplace transforms, respectively. 
In terms of the Fourier transform the Riesz-Feller fractional derivative in 
space is defined as 

J^{,D^f{xy,K} = -V,"(k) /(k) , = |Kre^(sign'^)^vr/2 , (2.4) 

where /(k) = J- {f{x); k} = f^^e~^^'^^f{x)dx. In other words is 
a pseudo-differential operatoiH with symbol xDq{k) = —iPq{k) , which is 
the logarithm of the characteristic function of the generic strictly stable (in 
the Levy sense) probability density, according to the Feller parameterization 



^ Let us recall that a generic pseudo-differential operator A, acting with respect to the 
variable a; € R, is defined through its Fourier representation, namely A /(k) — A{k) /(«:) 
where A{k,) is referred to as the symbol of A. The n-th derivative operator xD" = ^^jtt 
is a special case with symbol i,Z)"(k) — {—in)"^ . Generally speaking, a pseudo-differential 
operator A turns out to be defined through a kernel of a space convolution integral; 
this kernel is thus a sufficiently well-behaved function (absolutely) integrable in R which 
degenerates to a delta- type distribution when A — xD" . Furthermore, as a matter of fact, 
the symbol is given by the rule A{k) = (^e "*'*'") e+"°^ . 
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Figure 1: The Feller- Takayasu diamond of the generic Riesz- Feller derivative 



We note that the allowed region for the parameters a and 9 turns 
out to be a diamond in the plane {a , 6} with vertices in the points 
(0, 0) , (1, 1) , (2, 0) , (1, —1) , that we call the Feller- Takayasu diamona^ see 
Fig. 1. 

For a = 2 (henceforth = 0) we have xDq{k) = — = (—in)"^ , so we 
recover the standard second derivative. More generally for ^ = we have 

JD^{k) = = -(k2)"/2 so 




^0 = - -T^ • (2-5) 



In this case we refer to the LHS of (2.5) as simply to the Riesz fractional 
derivative of order a . We refer to the Appendix A for the explicit expression 
of the generic Riesz-Feller derivative. 



^Our notation for the stable distributions has been adapted from the original one by 
Feller. From 1998, see [23] . we have found it as the most convenient among the others 
available in the literature, see e.g. Janicki & Weron [31], Levy [33) . MontroU and associates 
[351 lii] . Samorodnitsky & Taqqu [H], Sato [5S], Uchaikin & Zolotarev [gSj, Zolotarev 
[67) . Furthermore, this notation has the advantage that all the class of the strictly stable 
densities are represented. As far as we know, the diamond representation in the plane 
{a, 9} was formerly given by Takayasu in his 1990 book on Fractals [59] . 
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Let us now consider the Caputo fractional derivative in time. Following 
the original idea by Caputo [3l H], see also [71 [8], [22], [50], a proper time 
fractional derivative of order (5 G (0, 1) , useful for physical applications, may 
be defined in terms of the following rule for the Laplace transform 

c[tDU{t);s} = s^Ks)-sP-^f{Q+), 0</3<l, (2.6) 

where f{s) = £{/(t);s} = e~^^ f{t) dt . Then the Caputo fractional 
derivative of f{t) turns out to be defined as 

{1 /■* f(i)(r) 
T{l-(3)Jo {t-rY (2.7) 

In other words the operator is required to generalize the well-known 
rule for the Laplace transform of the first derivative of a given (causal) 
function keeping the standard initial value of the function itself. We refer 
to the Appendix B for the relation of this derivative to the more common 
Riemann Liouville fractional integral and derivative. Here we report the 
most relevant formula (B.6) which provides alternative expressions for the 
Caputo fractional derivative for < /? < 1 : 



1 d /(r) - /(0+) 



r(i - /3) dt Jo {t-T)P ' r(i - 13) • 

It is worth to note that the time fractional derivative in the L.H.S. of Eq. 
(2.2) can be removed by a suitable fractional integration, see (B.4), leading 
to the alternative form 

1 /■* dr 

u{x, t) = u{x, 0+) + ^ xD^ u{x, r) _ ■ (2.9) 

Differentiating with respect to time we have another equivalent form 

It is well known that in the case of standard diffusion Eq (2.1) can be derived 
from the continuity equation, 

^^u{x,t) + ^F[u{x,t)]=0, (2.11) 
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where F is the flux given by 

F[u{x,t)] = -^u{x,t). (2.12) 

Whereas Eq. (2.11) is related to a conservation law of a physical quantity 
(therefore it has a universal character), Eq. (2.12) is a phenomenological law 
which states that the flux is proportional to the gradient of the field variable 
(the transported quantity) with opposite sign. It is met in several physical 
contexts with different names: e.g. when u is the temperature it is known 
as Fourier's law, when u is the pore pressure as Darcy's law, when u is a 
concentration of particles, as Pick's law. Here we use the last terminology in 
view of the possible applications in chemical physics. We recall that Pick's 
law is essentially an empirical law, which represents the simplest (local in 
space and time) relationship between the fiux F and the gradient of the 
concentration u observable in several physical phenomena. As a matter of 
fact, for some experimental evidences in complex transport phenomena, this 
law can be replaced by a more suitable phenomenological relationship which 
may account for possible non-local and memory effects, without violating the 
conservation law expressed by the continuity equation (2.11). Now it is not 
difficult to derive our fractional diffusion equation in the form (2.10) from a 
generalized Pick's law in which a suitable space-time operator depending on 
a, 9 and (3 is acting on the gradient. After simple manipulations based on 
recent results by Paradisi et al [48^ [l9] and Gorenfio et al [26] , we can write 



F{x,t) = tD'-^ ^Pi 



--u{x,t) 



(2.13) 



where tD^^^ denotes the Riemann-Liouville fractional derivative of order 
1 — /? (in time), see (B.l), and xPq is the pseudo-differential operator with 
symbol 

:= 4££M = |^|a-2gi(signK)e7r/2_ (2.14) 

Por a = 2 and /3 = 1 we recover the standard Pick's law since in this case 
jZ^o ^.p2 ^ J (Identity). Por 1 < a < 2 Eq. (2.13) results to be a a non- 
local connection between the fiux and the concentration gradient both from 
temporal and spatial view points. Por the nature of the operators involved, 
Eq. (2.13) can be referred to as the fractional Pick's law: it turns out to 
be consistent with the space-time fractional diffusion equation (2.2). We 
note that for < a < 1 Eq. (2.13) is meaningless since the symbol of the 
pseudo-differential Pq in (2.14) exhibits at k = a singularity not Pourier 
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integrable, which means that the kernel is not integrable in r[j. Generahzed 
Fick's laws with fractional derivatives have also been considered by other 
authors, including e.g. Zanette [65] and Caputo [5l|6]. 

3 The Green function for the space-time fractional 
diffusion 

The fundamental solution (or the Green function) for the space-time 
fractional diffusion is intended to be the solution of the governing equation 
(2.2), or (2.9) or (2.10), corresponding to the initial condition u{x^Q'^) = 
5{x) , It will be denoted by G^^(x,i). In the case of standard diffusion, 
see Eq. (2.1), the Green function is nothing but the Gaussian probability 
density function with variance a"^ = 2t , namely 

G°,i(x,t) = -^t-V2e-^V(4t). (3.1) 

In the general case, following the arguments by Mainardi, Luchko & Pagnini 
[37] we can prove that ^(j;, t) is still a probability density evolving in time 
with the noteworthy scaling property 

G%{x,t) = t-PI^K%(x/tPI'^) . (3.2) 

Here x/t^/" acts as the similarity variable and ^(•) as the reduced Green 
function. For the analytical and computational determination of the reduced 
Green function, from now on we restrict our attention to j; > because of 
the symmetry relation ^{—x) = K~^^{x) . Mainardi, Luchko & Pagnini 

^ From a purely mathematical view point one could overcome the above trouble for 
< a < 1 by stating the relationship between the flux F and the concentration u as 

F{x,t) = tD^-^ Hx,t)] , 

where xQe is the pseudo-differential operator with symbol 

:q|(^) := - 3m 3 W = IacI " - 1 e'(^^Sn^)ie + 1)^/2 

Physical reasons, however, lead us to avoid the range < a < f . In fact, for this 
range, xQ^i^-) would be a decreasing (or constant) function of |k| , which means that 
the contribution to the flux of the larger scales would be greater than (or equal to) that 
of the smaller ones, which is meaningless. 
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[37j have provided (for x > 0) the Melhn-Barnes integral representatiorj^ 



ax2mJ^-ioc T{1- ^s)T{ps)T{l- ps) ' 2a ' 

(3.3) 

where < 7 < minja, 1} . 

We recognize from the footnote that Eq. (3.3) by changing s into —s 
can be interpreted as a Melhn transform pair that allows us to write the 
Mellin transform of x ^{x) as 



a r(l + /3s/a) T{-ps) T{1 + ps)' 

(3.4) 

-min{a, 1} < K(s) < . 

In order to include s = in the convergence strip (so, in particular, the 
integral of K^ p{x) in Rq can be evaluated) we properly use in (3.4) the 
functional equation r(l -\- z) = zT{z) to obtain 

, ( ^ s. _ T{l-s/a)T{l + s/a)T{l + s) 

Jo """'/^^"^^^ ''''-''r{i-ps)r{i + ps)ni + ps/a)' 

(3.5) 

— min{a, 1} < 3ft(s) < a . 

In particular we find J^^°° ^(x) dx = p (with p = 1/2 if = 0). 

We note that Eq. (3.5) is strictly valid as soon as cancellations in the 
"gamma fraction" at the RHS are not possible. Then this equation allows 
us to evaluate (in Rq) the (absolute) moments of order 6 for the Green 



® The names refer to the two authors, who in the first 1910's developed the theory of 
these integrals using them for a complete integration of the hypergeometric differential 
equation. However, as pointed out in [11] (Vol. 1, Ch. 1, §1.19, p. 49), these integrals 
were first used by S. Pincherle in 1888. For a revisited analysis of the pioneering work of 
Pincherle (Professor of Mathematics at the University of Bologna from 1880 to 1928) we 
refer to the recent paper by Mainardi and Pagnini [38j . As a matter of fact this type of 
integral turns out to be useful in inverting the Mellin transform, as shown hereafter. If 

-M{/W;s} = r(s)= / f{r)r'-'dr, 7i<»(s)<72 



denotes the Mellin transform of a sufficiently well-behaved function /(r) , the inversion is 
provided by 



/'7-t-ioo 

M-'{r{s);r}^f{r) = — r{s)r-'ds 

J 7 — zoo 



where r > , 7 = 5R (s) , 71 < 7 < 72 ■ 
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function such that — minjo!, 1} < 6 < a . In other words, it states that 
p{x) = O as X — > +00. When cancellations occur in the 

"gamma fraction" the range of 6 may change. An interesting case is 
{a = 2, 9 = 0, < (3 < 1} (time-fractional diffusion including standard 
diffusion), where Eq. (3.5) reduces to 

+00 1 ■pi'i _|_ „\ 

„ At,W 2 f^lV^. ms)>-l. (3.0) 

This result proves the existence of all moments of order 6 > —1 for the 
corresponding Green function. In virtue of (3.2), (3.6) we have 

(3.7) 

and, for 5 = 2, the following formula for the variance 

/+°° n 2 a 

^x^G%{x,t)dx= ^^^^^^ tf^ , 0</3<l. (3.8) 

We note that the Mellin-Barnes integral representation allows us to construct 
computationally the fundamental solutions of Eq. (2.1) for any triplet 
{a, (3,9} by matching their convergent and asymptotic expansions, see |37j, 
|47j . Readers acquainted with Fox H functions can recognize in (3.3) the 
representation of a certain function of this class, see e.g. [29], [39] , [10] , [53] . 
[56] , [57] [58] , [62] . Unfortunately, as far as we know, computing routines 
for this general class of special functions are not yet available. 

Let us now point out the main characteristics of the peculiar cases 
of strictly space fractional diffusion, strictly time fractional diffusion, and 
neutral fractional diffusion based on the results stated in |37j . 

For /3 = 1 and < a < 2 {strictly space fractional diffusion) we recover the 
class of the strictly stable (non-Gaussian) densities exhibiting heavy tails 
(with the algebraic decay oc |x|~''""*'^'*) and infinite variance. For a = 2 
and < /? < 1 (strictly time fractional diffusion) we recover the class of 
the Wright-type densities exhibiting stretched exponential tails and finite 
variance proportional to t^ . Mathematical details on these two classes of 
probability densities can be found in [37j; for further reading we refer to 
Schneider [56] for stable densities, and to Gorenflo, Luchko & Mainardi 
[201 [5T] for the Wright- type densities. 

As for the stochastic processes governed by these distributions we can expect 
the following. 
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For the case of non-Gaussian stable densities we expect a special class 
of Markovian processes, called stable Levy motions, which exhibit infinite 
variance associated to the possibility of arbitrarily large jumps {Levy flights), 
whereas for the case of Wright-type densities we expect a class of stochastic 
non-Markovian processes, which exhibit a (finite) variance consistent with 
slow anomalous diffusion. 

For the special case a = P < 1 (neutral difl^usion) we obtain from (3.3) an 
elementary (non-negative) expression 



where < 7 < a . 

For the generic case of strictly space-time diffusion (0 < a < 2, < < 1), 
including neutral diffusion for a = P < 1, Mainardi, Luchko & Pagnini 
[37j have proven the non negativity of the corresponding reduced Green 
function and consequently its interpretation as probability density. In this 
case we obtain a class of probability densities (symmetric or non-symmetric 
according to ^ = or ^ 7^ 0) which exhibit heavy tails with an algebraic 
decay oc . Thus they belong to the domain of attraction of the 

Levy stable densities of index a and can be referred to as fractional stable 
densities, according to a terminology proposed by Uchaikin [60]. The 
related stochastic processes are expected to possess the characteristics of 
the previous two classes. Indeed, they are non-Markovian (being /3 < 1) 
and exhibit infinite variance associated to the possibility of arbitrarily large 
jumps (being a < 2). 

4 The discrete random walk models for the 
Markovian fractional diffusion 

It is known that a numerical approach to the standard diffusion equation 
(2.1) based on a proper finite difference scheme provides a discrete 
Markovian random walk model for the classical Brownian motion, see e.g. 
Zauderer j66j . 

In this section we intend to generalize this approach (that will be hereafter 
recalled) in order to provide a discrete Markovian random walk model for 




(3.9) 




X > 
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the Levy stable motion of any order a G (0, 2) and skewness 9 restricted as in 
(2.3). For this purpose we present a notable finite-difference approach to the 
strictly space fractional diffusion equation subjected to relevant restrictions, 
as we shall show in the following. 

The common starting point of our analysis is obviously the discretization of 
the space-time domain by grid points and time instants as follows 

Xj=jh, h>0, j = 0,±1,±2,... ; ,^ ^, 

tn = nT, T>0, n = 0,l,2,... 

where the steps h and r are assumed to be small enough. The dependent 
variable u is then discretized by introducing yj{tn) as 

rXj+h/2 

Ujitn)^ I u{x,tn)dx K, hu{Xj,tn) ■ (4.2) 

Jxj~h/2 

4.1 The discrete random walk model for the standard 
diffusion 

Let us now consider the standard diffusion equation (2.1). With the 
quantities yj{tn) so intended, we replace Eq. (2.1), after multiplication by 
the spatial mesh-width h, by the finite-difference equation 

yj{tn+l) - yjjtn) ^ yj+l{tn) - '^yjjtn) + yj-l{tn) . . o^ 

accepting that for positive n in (4.3) we have approximate instead of exact 
equality. Since we are interested to approximate the fundamental solution 
(the Green function), we must equip (4.3) with the initial condition yj(0) = 
5j , where the Kronecker symbol represents the discrete counterpart of the 
Dirac delta function. 

This approach can be interpreted as a discrete (in space and time) 
redistribution process of some extensive quantity provided it is conservative. 

If the extensive quantity is non-negative, e.g. mass or a sojourn probability, 
we have to preserve its non-negativity. In the first case the yj{tn) are 
imagined as clumps of mass, sitting at grid points x = Xj in instants t = tn, 
which collect approximatively the total mass in the interval Xj — h/2 < x < 
Xj + h/2 . In the second case, the yj{tn) may be interpreted as the probability 
of sojourn in point Xj at time tn for a particle making a random walk on 
the spatial grid in discrete instants. From now on, we agree to pursue the 
probabilistic point of view. 
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In order to have a conservative and non-negativity preserving redistribution 
process, the discrete variable yj is subjected to the conditions 

+00 +00 

E yji^n) = E > yji^n) > , for all i G Z , n G Nq . (4.4) 

j=—oo j=—oo 

We easily recognize that our discrete redistribution process is akin to a 
Markov chain: when time proceeds from t = t„ to t = tn+i , the sojourn- 
probabilities are redistributed according to the transition law 

00 

Vjitn+i) = Pk yj-k{tn) , J e Z , n G No , (4.5) 

k=—oo 

where the denote suitable transfer coefficients, which represent the 
probability of transition from xj^k to xj (likewise from Xj to Xj^k)- The 
transfer coefficients are to be found consistently with the finite-difference 
equation (4.3) equipped with the proper initial condition. The process turns 
out to be both spatially homogeneous (the probability of jumping from a 
point Xj to a point xj^^ ^ot depending on j) and time-stationary (the pk not 
depending on n), as is advised when considering our Cauchy problem and the 
definition of the difference operators. Furthermore the transfer coefficients 
must satisfy the conditions 

oo 

Pk = l, Pk>0, k = 0,±l,±2, ... (4.6) 

A;=— oo 



The transfer coefficients in our special case are easily deduced from (4.3) 
and (4.5): they turn out to be 

Po = l-2^, p±i = ^, p±k = 0, fc = 2,3...,. (4.7) 

subject to the condition 

r 1 

We refer to Eqs (4.6)-(4.8) as the basic equations for the standard random 
walk model for the Gaussian process. The constant /i will be denoted as 
the scaling parameter of the standard diffusion equation. We recognize that 
only jumps of one step to the right or to the left or jumps of width zero 
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Figure 2: Histogram (left) and a sample path with increments (right) for 
standard diffusion {a = 2, /3=1, 9 = 0} (Brownian motion) 

occur. This corresponds to a well-known simple approximate realization of 

the Brownian motion. 

The finite difference scheme (4.3) can be used for producing sample paths 
of individual particles performing the Brownian motion and for producing 
histograms of the approximate realization of the Gaussian density, by 
simulating many individual paths with the same number of time steps and 
making statistics of the final positions of the particles. 

Our simulations, based on 10,000 realizations, have been limited to the 
interval |x| < 5 , where we have considered the random walks to take place, 
have produced an histogram for the Gaussian density at t = 1 , sec Fig 2 
(left plate). Then, particles leaving this space have been ignored. In Fig 2 
(right plate) we have displayed a particular sample path or trajectory (up) 
and the corresponding increment series (below) obtained with N = 500 time 
steps, by taking an intermediate value for the scaling factor (/x = 0.4 < 1/2) 
and reasonable values for the space and time steps (/i = 0.0707 , r = 0.002) 
in order to get Nt = t = 1 . Of course the histogram has been depicted 
by adopting a space step much larger than for the sample path, namely 
h = 0.25 . 
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4.2 A discrete random walk model for the strictly space 
fractional diffusion 



Discretizing all the variables as done for the standard diffusion equation, 
namely introducing a space-time mesh of widths h and r and a discrete 
variable yj{tn) interpreted as in (4.2), the essential idea is to replace the 
strictly space fractional diffusion equation by the finite-difference equation 



where the difference operator hDg is intended to converge to Dq (discussed 
in Appendix A) as /i — > . As usual, we have adopted a forward difference 
quotient in time at level t = tn for approximating the first-order time 
derivative. For ^Dq we require a scheme which must reduce as a = 2 
to a symmetric second-order difference quotient in space at level t = 
tn , which has been adopted for approximating the second-order space 
derivative. Furthermore, the the finite-difference equation (4.9) is required 
to be consistent with a conservative, non-negativity preserving redistribution 
process, subject to the initial condition yj(0) = 5jo . Referring to Eqs (A. 4), 
(A. 8) and (A.lO)-(A.ll), we write 



With the notation hD± we intend to adopt hD-\- when < ^ < 1 and hD- 
when -1 < ^ < . 

We thus must find suitable finite-difference schemes for the pseudo- 
differential operators entering the Feller derivative, namely hD± (Weyl 
fractional derivatives of order a) if a 7^ 1 and /i-Dg (Riesz derivative of 
first-order) if a = 1 . For this purpose we must keep distinct the three cases 



For a 7^ 1 the starting point is the Griinwald-Letnikov discretization of 
fractional derivatives, on which the reader can inform himself from the 
treatises on fractional calculus, see e.g. [32], [56], [50], [53], or in the review 
article by Gorenflo [18]. However, for our purposes, we must make a clever 
use of the Griinwald-Letnikov scheme treating separately the two cases 



hD^yjitn), < a < 2,|0| < min{Q,2-Q}, (4.9) 



r 




(4.10) 




(4.11) 
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(a) and (b) in order to be consistent with a conservative, non-negativity 
preserving redistribution process. We have 



1 oo I, /q;\ 

^ E (-1) L 1 ViTk , in the case (a) < a < 1 , 

I oo I, /a\ 

[j.) ^iiiTfc ' in the case (b) 1< a < 2 . 



(4.12) 

Notice the shift of index in the case (b) which among other things has the 
effect that in the hmiting case a = 2 (the classical diffusion equation) we 
obtain the standard symmetric three-point difference scheme. For more 
details and discussions see Gorenflo & Mainardi |23j . 

For a = 1 we limit ourselves to consider the case 9 = 0, that requires 
the discretization of the Hilbert transform in (A. 11). The hyper-singular 
integral representation of the symmetric space-fractional derivative (the 
"Riesz derivative") given by Gorenflo & Mainardi in |25| [as formula (2.20)] 
readily offers us the discretization 

'^^y^ - Vh — w^) — ■ ^ ^ ^ 

For < 1^1 < 1 we have analogously to (A. 10) to take into account the 
forward/backward difference quotients for hD± 

hD^-^^, hD.-^^, (4.14) 

to obtain a non-negative scheme in a proper way. However, we here leave 
out these "drifts" (in negative or positive directions) from the computation 
of the transfer coefficients pt- 

Then, inserting the expressions (4.12)-(4.13) in the finite-difference equation 
(4.9)- (4. 10) yields the transition law 

oo 

yj{tn+i) = XI PkVj-kitn), jeZ, nGNo, (4.15) 

fc=— oo 

where the transfer coefficients pk turn out to be, respectively, 

^ cos (6'7r/2) 

^ cos (011/2) ' , , 

^ ' ^ < a < 1, |6l| < a; 




c±, k = l,2, 



(4.16a) 
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Pq = 1 + fi 

p±l = -At 
p±k = (-I)V 
1 



(c+ + c_) = l — /ia 



cos(6l7r/2) 
I cos (a7r/2)| ' 

l<a<2, 16*1 <2-a; 



(4.166) 



a 
+ 



2,3,, 



Po 

P±k = - 



2/u 



1 



a 



1, 61 = 0; 



vr k{k + 1) ' 
with the scale parameter 



1,2,... 



n := — , < a <2, 



(4.16c) 



(4.17) 



It is straightforward to check in all cases the summation condition 

oo 

Pk = \ . Since all p±k ^ , A; 7^ the non-negativity condition is 

k=—oo 

met if we require pq > 0, i.e. if the scale parameter /i is restricted as follows 

' cosQ7r/2 1 i/ii / 

< a < 1, ^ < a, 

cos t^7r/2 

1 I cos a7r/2| , ^, /a-i q\ 
^4r' l<a<2,^<2-Q, (4-18) 

-, a = i, e = o. 

L 2' 



T 

l-l ■ — ^ fJ'max 



We note that in both the limits a — > 1~ and a — > 1"*" the permissible range 
of the scaling factor fi is vanishing. In numerical practice the consequence 
will be that if a is near 1 the convergence is slow: for good approximation 
we will need a very small step-time r with respect to the step-length 



The striking difference to our discretized Brownian is the appearance of 
arbitrarily large jumps with power- like decay probability for which these 
discrete models are referred to as Levy flights. See Fig. 3 for a sketch of the 
transition scheme and Section 6 for the numerical results. 

^In order to get a continuous transition to the case q = 1 we need to consider a different 
discretization scheme, presented in [T2] for the first time and rigorously analyzed in [25] . 
For this scheme, however, we loose the continuity as a — > 2~ . This means that "there is 
no free lunch"; we have to pay for the good behaviour at a = 1 with bad behaviour at 
a = 2 in a sense described in [19j . 
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I ' ' ' ' ' ' ' ' ' 1 

-j ... -2 -1 O 1 2 ... j 

Figure 3: Sketch of the transition scheme for the space- fractional random 
walker 

5 The discrete random walk models for the non- 
Markovian fractional diffusion 

Let us now try to generalize the above arguments by adapting them to the 
non-Markovian cases of our space-time fractional diffusion equation (2.2). 
We find it convenient to proceed by steps: we first consider the case of 
the strictly time fractional dijfusion, see also [26j, and then, by combining 
these arguments with those for strictly space fractional diffusion, we treat 
the strictly space-time fractional dijfusion. 

5.1 A discrete random walk model for the strictly time 
fractional diffusion 

Discretizing all the variables as for the standard diffusion equation, namely 
introducing a space-time mesh of widths h and r and a discrete variable 
yj{tn) interpreted as in (4.2), the essential idea is to replace the strictly time 
fractional diffusion equation by the finite-difference equation 

.Z^^.(Wi) = ^^-"^^^"^"^^;i^-^+"^--^^^"\ </?<!, (5.1) 
where the difference operator t-D^ is intended to converge to as r — > . 
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As usual, we have adopted a symmetric second-order difference quotient in 
space at level t = tn for approximating the second-order space derivative. 
For T-D^ we require a scheme which must reduce as /? = 1 to a forward 
difference quotient in time at level t = tn , which is usually adopted for 
approximating the first-order time derivative. Then, for approximating 
the time fractional derivative (in Caputo's sense), wc adopt a backward 
Griinwald-Letnikov scheme in time (starting at level t = tn+i) which reads 



ra+l 



fc=0 



< /3 < 1 . (5.2) 



Here the subtraction of tjj (0) in each term of the sum reflects the subtraction 
of /(O"*") in formula (B.6) for the Caputo fractional derivative. Combining 
(5.1) and (5.2), introducing the scaling parameter 



/x:=-p-, </?<!, 



(5.3) 



and using the "empty sum" convention ^ = if g < p (here p = 1 when 

k=p 

q = n = 0), we obtain for n > (to = 0): 



\k+l 



yj{tn+l-k) 



k=0 



k=l 



(5.4) 



Thus, (5.4) provides the universal transition law from tn to tn+i valid for all 
n > 0. For convenience let us introduce the coefficients Cfe, bm 



Ck = (-1) 



k+i 



k>l, 



E(-i) 



fe=0 



, k j 



(5.5) 



m > 0. 



For /3 = 1 (standard diffusion) we note that all these coefficients are 
vanishing except 6o = ci = 1 . For < /? < 1 they possess the properties 



^ Cfc = 1 , 1 > /3 = Ci > C2 > C3 > . . . 

fe=i 



0, 



(5.6) 
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(5.7) 



bo = l = ^ Ck, 5^ = 1 - ^ Cfc = J2 ^k, 

k=l k=l k=m+l 

I = bo > bi > b2 > > . . . ^ . 

We thus observe that the and the bm form sequences of positive numbers, 
not greater than 1, decreasing strictly monotonically to zero. Thanks to the 
introduction of the above coefficients the universal transition law (5.4) can 
be written in the following noteworthy form 

n 

yj{tn+l) = bn Vjito) + ^Ck yj{tn+l-k) + M [?/j+l(*n) - '^Vjitn) + %-l(*n)] , 
k=l 

(5.8) 

with the empty sum convention convention if n = . In particular we get, 
for n = : 

%(*i) = (1 - Vjito) + 1^ [yj+i{to) + yj-i(to)] ; 

for n = 1 : 

yj{t2) = biyj{to) + (ci - 2//) yj{ti) + fi [yj+i{h) + yj-iih)] ; 
forn > 2 : 

n 

Vjitn+l) = bn yjito) + Cfc yj{tn+l-k) 

k=2 

+(ci - 2/x) yj{tn) + H [yj+i{tn) + %-l(in)] • 

Observe that C\ = j3 . The scheme (5.8) preserves non-negativity, if all 
coefficients are non-negative, hence if 



Furthermore it is conservative, as we shall prove by induction: i.e. 



(5.9) 



+00 +00 +00 

E \yjito)\ < oo ^ E = E %(*o) > " e N. (5.10) 

j=—oo j=—oo j=—oo 

+ 00 

In fact, putting Sn = E Vjitn) for n > , then from (5.8) we get 

j=-oo 

5i = (1 - 2fi) ^yj(io) + /xE%-i(*o) +y"E%+i(*o) = ^o,, 
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and for n > 1 we find always from (5.8), assuming Sq = Si = . . . Sn already 
proved, 

n 

Sn+l =bnSo + ^Ck Sn+l-k + {P - 2fl + fl + fi) Sn 
k=2 

= (^K + cfc^ So = So, 

using P = ci . We have thus proved conservativity. Non-negativity 
preservation and conservativity mean that our scheme can be interpreted as 
a redistribution scheme of clumps yj{tn) ■ For orientation on such aspects and 
for examples let us quote the works by Gorenflo on conservative difference 
schemes for diffusion problems, see e.g. [T6l [T7] . 

The interpretation of our redistribution scheme is as follows: the clump 
yj{tn+i) arises as a weighted-memory average of the (previous) n + 1 values 
Ujitm) , with m = n , n — 1 , . . . ,1,0, with positive weights 

n 

(3 = Ci, C2 , . . . ,Cn , 6„ = l-^Cfe, (5.11) 

k=l 

followed by subtraction of 2^yj{tn), which is given in equal parts to 
the neighbouring points and Xj+i but replaced by the contribution 

/i[yj+i(in) + Vj-iitn)] from these neighbouring points. For random walk 
interpretation we consider the yj{tn) as probabilities of sojourn at point Xj 

+ 00 

in instant t„ requiring the normalization condition ^ Vjito) = 1 • 

j=-oo 

For n = Equation (5.8) means (by appropriate re-interpretation of the 
spatial index j): A particle sitting at xj in instant to jumps, when t 
proceeds from to to ti , with probability fi to the neighbour point xj+i , 
with probability fi to the neighbour point xj-i , and with probability 1 — 2// 
it remains at Xj . For n > 1 we write (5.8), using /3 = ci , as follows: 

yj{tn+l) = (^^-Yl '^'^^ + ^nVjih) + C„_iyj(t2) + . . . + C2yj{tn-l) 

+{ci - 2fl) yj{tn) + fJ, [yj+l{tn) + yj-l{tn)] . 

(5.12) 

Obviously, all coefficients (probabilities) are non negative, and their 
sum is 1. But what does it mean? Having a particle, sitting in xj 
at instant tn , where will we find it with which probability at instant 
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t„_l_i? Prom (5.12) we conclude, by re-interpretation of the spatial index 

j, considering the whole history of the particle, i.e. the particle path 
{x(to) , x{ti) , x{t2) , • • • , x{tn)} , that if at instant t„ it is in point xj , there 
is the contribution ci — 2// to be again at xj at instant tn+i , the contribution 
fx to go to Xj-i , the contribution /x to go to xj^i . But the sum of these 
contributions is ci = /? < 1 . So, excluding the case /? = 1 in which we 
recover the standard diffusion (Markovian process), for (3 < 1 we have to 
consider the previous time levels {non- Markovian process). Then, from level 
tn-i we get the contribution C2 for the probability of staying in xj also at 
time tn+i , from level tn-2 we get the contribution C3 for the probability of 
staying in xj at time tn+i , ■ ■ ■ , from level ti wc get the contribution c„ for 
the probability of staying in Xj at time tn+i , and finally, from level to = 
we get the contribution 6„ for the probability of staying in Xj at time tn+i ■ 
Thus, the whole history up to tn decides probabilistically where the particle 
will be at instant tn+i ■ 

Let us consider the problem of simulation of transition from time level tn 
to tn+i- Assume the particle sitting in Xj at instant i„ . Generate a random 
number equidistributed in < p < 1 , and subdivide the interval [0, 1) as 
follows. Prom left to right beginning at zero we put adjacent intervals of 
length ci,C2, . . . ,c„,6„, for consistency left-closed, right-open. The sum of 
these is 1. We divide further the first interval (of length ci) into sub-intervals 
of length /X, ci — 2/x, /x . Then we look into which of the above intervals the 
random number falls. If in first interval with length ci = fi + (ci — 2^) + (i, 
then look in which subinterval, and correspondingly move the particle to 
Xj-i , or leave it at xj or move to Xj+i . If the random number falls into one 
of the intervals with length C2 , C3 , . . . , c„ (ie. with 2 < k < n), then 
move the particle back to its previous position x{tn+i-k), which by chance 
could be identical with xj = x{tn) ■ If the random number falls into the 
rightmost interval with length then move the particle back to its initial 
position x{tQ) , for which we recommend x{tQ) = 0, meaning yj{to) = Sjo , 
in accordance with the initial condition u{x,0) = S{x) for (2.2) with a = 

2,e = o. 

A sketch of the transition scheme for the random walker is reported in Figure 
4. Besides the diffusive part {/j, , ci — 2/x , fj,) which lets the particle jump at 
most to neighbouring points, we have for < /3 < 1 the memory part which 
gives a tendency to return to former positions even if they are far away. Due 
to Equations (5.6)-(5.7), of course, the probability to return to a far away 
point gets smaller and smaller the larger the time lapse is from the instant 
when the particle was there. 
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Figure 4: Sketch of the transition scheme for the time-fractional random 
walker 

5.2 A discrete random walk model for the strictly space-time 
fractional diffusion 

By combining the approach of the preceding Subsections 3.2 and 4.1 we 
can construct a discrete random walk model for the strictly space-time 
fractional diffusion equation, namely for the casejO <q;<2,0</3<1}. 
We replace in (2.2) the (Riesz-Feller) space-fractional derivative by (4.10) 
and the (Caputo) time-fractional derivative by (5.2). Solving then for the 
"new" value yj{tn+i) we see that the scaling parameter ^ = t^^ /h'^ plays 
the essential role in obtaining a scheme with all transition coefficients non- 
negative. In fact, also here jj, must be restricted to a suitable interval 
(0, Umax] ■ In analogy to the case of the strictly time fractional diffusion 
the "discrete" diffusion in space occurs only between the time-level t„ and 
tn+i and the memory part of the process only straight-backwards in time. 
However, in contrast to the strictly time fractional diffusion, the discrete 
diffusion (or the random walker) can now go to any grid point in space, not 
only to immediate neighbouring grid points. We abstain from presenting 
here all the lengthy formulas that again need distinction between the three 
cases: (a) < a < 1 , (b) 1 < a < 2 , and (c) a = 1 . 
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6 Numerical results and concluding discussions 



Our simulations are all based on 10,000 realizations. In addition to the 
Brownian motion, see Fig. 2, we have considered a set of 8 case-studies, 
see Figs 5-12, for our space-time fractional diffusion that, in our opinion, 
can better illustrate the state-of art of this analysis. The sample paths and 
the corresponding increments are plotted against the time steps up to 500, 
see the right figure-plates, while the histograms refer to densities at t = 1 
for \x\ < 5 , see the left figure-plates. All the plots were drown by using 
the MATLAB system. The relevant parameters a, 6, f3, n, h and r used in 
the Figures are reported in Table I. For convenience we have also reported 
the maximum value of the scaling parameter fi as can be deduced from our 
theoretical analysis. More details can be found in [45j . 
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Table I: The relevant parameters for the simulations 
a = space-fractional order, 9 = skewness, (3 = time-fractional order, 
jjL = T^//i" = scaling parameter, 
^5 = space-step, ts = time-step for sample paths, 
hn = space-step, th = time-step for histograms. 



In practice, in our numerical studies there is required truncation for two 
different causes. As in the classical Brownian motion, a trivial truncation 
is required if a priori one wants a definite region of space to be considered 
in which the walk takes place. Then, particles leaving this space have been 
ignored. However, if < a < 2 , at variance with our discretized Brownian 
motion, we now have an infinite number of transition probabilities. Since 
it is impossible to simulate all infinitely many discrete probabilities, so the 
size of possible jumps must be limited to a maximal possible jump length. 
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The plates in Fig. 5 are concerning two cases of symmetrical, strictly space 
fractional diffusion: {a = 1.75, P = 1, 9 = 0}, {a = 1.5 , /3 = 1, 6* = 0} 
whreas the plates in Fig. 6 are concerning two cases of extremal, strictly 
space fractional diffusion: {a = 1.75, (3 = 1, 9 = —0.25}, {a = 1.5, /? = 
1, 9 = —0.50} From the sample paths in Figs. 5, 6 one can recognize the 
"wild" character (with large jumps) of the Levy flights with respect to the 
"tame" character of the Brownian motion outlined in Fig. 2. 




I _Q 4 1 1 1 1 1 1 1 1 1 1 1 

-5 -4 -3 -2 -1 1 2 3 4 5 '0 100 200 300 400 500 



Figure 5: Histograms (left) and sample paths with increments (right) for 
symmetrical, strictly space fractional diffusion. 

The plates in Fig. 7 are concerning two cases of strictly time fractional 
diffusion: {a = 2, (3 = 0.75} and {a = 2, (5 = 0.50}; here the paths exhibit 
the memory effect visible in a kind of stickiness combined with occasional 
jumps to points previously occupied, in distinct contrast to the rather tame 
behaviour in case {a = 2, P = 1} (simulation of Brownian motion). 
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Finally, the plates in Fig. 8 are concerning two cases of strictly space-time 
fractional diffusion: {a = 1.50, /? = 0.50, 9 = 0} and {a = 1.50, (3 = 
0.50, 6 = —0.50}, where the combined effects of the previous cases are 
present. 

We have used our discrete models for simulation of particle trajectories by 
interpreting our redistribution schemes as descriptions of Markov chains 
with infinitely many states, namely the possible positions xj. However, as 
they have the form of difference schemes in a regular space-time grid they 
can "in principle" also serve for the purpose of approximate computation of 
the temporal evolution of the density u{x, t). Namely, we are expecting that 
u{xj,tn) is approximated by yj(tn)/h. In fact, from their essential properties 
(conservativity and preservation of non- negativity, see (4.6), (5.11) and 
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Figure 7: Histograms (left) and sample paths with increments (right) for 
strictly time fractional diffusion. 

(5.12) ) it can by standard methods of numerical analysis been shown 
that our models interpreted as difference schemes are stable and consistent, 
hence convergent. We say "in principle" because for practical application 
appropriate truncations are required. That there is also convergence in 
the chain interpretation (complete convergence in the stochastic sense) can 
be shown by methods of Fourier and Laplace analysis, see Gorenfio and 
Mainardi in [24:\ [25j for the space fractional diffusion. 

There are other methods of simulating space-fractional diffusion processes. 
Without attempting to be exhaustive, without going into details and without 
attempting to give a survey of the existing and growing literature on the 
subject let us mention a few possibilities: generations of random numbers 
distributed according to a given stable law, random walks discrete in time 
but continuous in space (namely jumps to arbitrary places in space, not only 
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Figure 8: Histograms (left) and sample paths with increments (right) for 
strictly space-time fractional diffusion. 

to grid points). Furthermore, let us mention here the Chechkin-Gonchar 
random walk (suggested and used in [9] and mathematically analyzed in 
[25j). and simulation via compound Poisson processes (with proper scaling 
of space and time), see [27] . 

In the special field of methods, discrete in space and time, there also are 
other methods available (see the methods sketched in [19] and discussed 
in detail in [241 ES] of which, in particular, the Gillis- Weiss random walk 
must be cited, see [H]). Final word: Still many questions are open in this 
challenging and fascinating field of research! 
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Appendix A: The Riesz-Feller space fractional derivatives 

In this Appendix we provide the explicit expression of the Riesz-Feller 
fractional derivative xDg which, according to (2.4) is defined as the pseudo- 
differential operator with symbol 

Let us now express more properly the operator — xDg as inverse 
of a suitable integral operator ^Iq , whose symbol is required to be 
e ~^ (^^Sn ^) ^^/2 ^ go we may write 

xDq := — • (^-2) 

This integral operator was found by Feller [T2] in 1952 generalizing the 
approach by Marcel Riesz to Fractional Calculus, see e.g. [51J, and it is 
referred to as Feller potential by Samko, Kilbas Sz Marichev [53j^ . Using 
our notation, the Feller potential reads 

fix) = c_(a, e) Jl /(^) + 0) f{x) , (^.3) 

where, if0<a<2, a^l, 

sin [{a - 9) 7r/2] sin [(a + 9) 7r/2] 

c+[a,9) = — — , c_(a,6') = — — r , (^.4) 

sm (avrj sm(a7rj 

and, by passing to the limit (with 6* = 0) c+(2,0) = c_(2,0) = -1/2. In 
(A. 4) the operators xl± denote the Riemann-Liouville fractional integrals, 
also known as Weyl fractional integrals) , defined as 

./?/(x) = -i- r(x-o"-V(e)de, 



1 r+oo 

■I1f{x) = -— / {i-x)^-^ f{£,)dC 
F a Jx 



\a) 

We note that whereas the coefhcients c± can loose their meaning when a is 
an integer, the Riemann-Liouville integral operators are well defined in 



We must note that in his original paper FeUer used a skewness parameter 5 different 
from our 9 ; the symbol of the Feller potential is 



(l^le-^ (sign n)5 



so — , 6' = ao . 

2 a TT 



Feller and Samko, Kilbas & Marichev thus use where their & is related to our Q as 
above. 
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their action on rapidly decreasing integrable functions, for any a > , being 
set equal to the identity operator when q = , for convenience (justified by 
passage to the limit a ^ 0). In the particular case 9 = we get 

c+(a,0)=c_(a,0) = - ^— - , (A6) 

2 cos [TTa/2) 

and thus we recover the Riesz potential, see e.g. [53] , 

JS /(^) := , \ - Cr-i /(O . iA.7) 

2T[a) cos(7ra/2j J-oo 

The Riesz and the Feller potentials are well defined if the index is located 

in the interval (0 , 1) , or in the interval (0, 2) , and we have the semigroup 

property, xlg xIq = xlg'^^ ,ifO<a<l, 0</3<l and a + (3 < 1 . Then, 

following Feller, we define by analytic continuation the pseudo-differential 

operator (A. 2) in the whole range 0<a<2,a7^1,as 

xD'^ ■.= -[c+{a,e) xD1 + c.{a,9) xD'^] , 0<a<2. (A8) 

where xD'^ ■= xl^"' and ^-D" := x^Z^ are the inverses of the operators a,/" 
and xl- ) respectively, and are referred to as the Weyl fractional derivatives. 
When = the Riesz-Feller derivative can be simply referred to as Riesz 
derivative. We note from (A. 4) the property c±{a,9) = c^(— q). For 
integral representations of the operators see [53]; we have 

[±-^ (Jl-") . it 0<£.< 1, 

For a = 2 (0 = 0) we recover the standard second derivative. In fact, 
from (A.6) c+(2,0) = c_(2,0) = -1/2, so from (A.8)-(A.9) xD^ = -I^^ = 
{xl^' + xlZ') /2={£ + £)/2 = i,= xD' . For a = 1 (1^1 < 1) we 
have (from calculation with symbols) 

xDl fix) = [cos(07r/2) xDl + sin(07r/2) f{x) , (A.W) 

where xDlf{x) = - xD [xH f{x)] , with 

vr J-oo X - t, TT J-oo c, 

In (A. 11) xH denotes the Hilbert transform (with symbol xH{k) = isign k ) 
and its singular integral is understood in the Cauchy principal value sense. 
We note that in the limiting extremal cases 9 = ±1 we recover the standard 
first derivative, i.e. = ± xD . 
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Appendix B: The Riemann-Liouville and Caputo time frac- 
tional derivatives 



The purpose of this Appendix is to clarify for the interested reader the 
main differences between the Caputo fractional derivative t-D* adopted in 
this paper, see (2.8), and the Riemann-Liouville fractional derivative 
usually adopted in the literature. Any formula here is valid for t > , with 
the assumption that the function f{t) has a finite limit as t — > 0+ . We recall 

The two fractional derivatives are related to the Riemann-Liouville fractional 
integral as follows. The Riemann-Liouville fractional integral is 

tJ" fit) ■■= ^ /J fir) it - r)^'-i dr, ^i > , {B.2) 

(with the convention f J*^ f{t) = fit)) and is known to satisfy the 
semigroup property tJ^tJ'^ = tJ^^'^ , with > 0. For any (3 > the 
Riemann-Liouville fractional derivative is defined as the left inverse of the 
corresponding fractional integral (like the derivative of any integer order), 
namely tL*^ * f{t) = f{t) . Then for /3 e (0, 1] we have 

tD" fit) ■.=t tJ'~f fit) , tD^^fit) :=t j'-f tD' fit) , {B.3) 

tjP tD^, fit) = tJ^ tJ'-P tD' fit) = tJ' tD' fit) = /(t) - /(0+) . iB.4) 

Recalling the rule 



(B,5) 



it turns out for < /? < 1 , 

fit) = tD^ [fit) - /(0+)] = tD^ fit) - /(0+) ^^^^ . (i?.6) 

Note that for [3 = 1 the two types of fractional derivative coincide. 
The Caputo fractional derivative represents a sort of regularization in the 
time origin for the Riemann-Liouville fractional derivative and satisfies the 
relevant property of being zero when applied to a constant. For more details 
on this fractional derivative (and its extension to higher orders) we refer the 
interested reader to Gorenflo and Mainardi [22], Podlubny 'SO] and Butzer 
and Westphalj2j. 
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